using System;
using System.Numerics;
using FFTWSharp;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;


namespace FFTWSharp
{
    public class FFTWWrapper
    {
        private fftwf_plan plan;
        private fftwf_complexarray mfin, mfout;

        public FFTWWrapper()
        {

        }

        /// <summary>
        /// 计算功率谱（RMS）
        /// </summary>
        /// <param name="input">输入信号</param>
        /// <param name="fs">采样频率</param>
        /// <param name="amp">幅值</param>
        /// <param name="pha">相位</param>
        /// <param name="coorX">横坐标</param>
        /// <returns></returns>
        public bool FFTPowerCompute(float[] input, float fs, ref double[] amp, ref double[] pha, ref double[] coorX)
        {
            fftwf.fftwf_forget_wisdom();

            InitialPlan(input);
            int N = input.Length;
            float df = 1;
            if (N != 0)
                df = fs / N;
            float[] hanning = Hanning(input);
            float[] i2 = new float[N * 2];
            for (int i = 0; i < N; i++)
            {
                i2[2 * i] = hanning[i];
                i2[2 * i + 1] = 0;
            }

            mfin.SetData(i2);
            plan.Execute();
            Complex[] complex = mfout.GetData_Complex();

            amp[0] = complex[0].Magnitude / N * (1.633 / Math.Sqrt(2));
            pha[0] = complex[0].Phase * 180 / Math.PI + 90;
            coorX[0] = 0;

            for (int i = 1; i < N / 2; i++)
            {
                amp[i] = complex[i].Magnitude * (2.0 / N) * (1.633 / Math.Sqrt(2));
                pha[i] = complex[i].Phase * 180 / Math.PI + 90;
                //double p = Math.Atan(complex[i].Imaginary / complex[i].Real);
                //pha[i] = p * 180 / Math.PI + 90;
                coorX[i] = df * i;
            }

            return true;
        }

        /// <summary>
        /// 计算幅值谱（RMS）
        /// </summary>
        /// <param name="input">输入数据数组</param>
        /// <param name="fs">采样频率</param>
        /// <param name="amp">幅值</param>
        /// <param name="pha">相位</param>
        /// <param name="coorX">横坐标</param>
        /// <returns></returns>
        public bool FFTAmpCompute(float[] input, float fs, ref double[] amp, ref double[] pha, ref double[] coorX)
        {
            fftwf.fftwf_forget_wisdom();

            InitialPlan(input);
            int N = input.Length;
            float df = 1;
            if (N != 0)
                df = fs / N;
            float[] hanning = Hanning(input);
            float[] i2 = new float[N * 2];
            for (int i = 0; i < N; i++)
            {
                i2[2 * i] = hanning[i];
                i2[2 * i + 1] = 0;
            }

            mfin.SetData(i2);
            plan.Execute();
            Complex[] complex = mfout.GetData_Complex();

            amp[0] = complex[0].Magnitude * (1.0 / N) * (2 / Math.Sqrt(2));
            pha[0] = complex[0].Phase * 180 / Math.PI + 90;
            coorX[0] = 0;

            for (int i = 1; i < N / 2; i++)
            {
                amp[i] = complex[i].Magnitude * (2.0 / N) * (2 / Math.Sqrt(2));
                pha[i] = complex[i].Phase * 180 / Math.PI + 90;
                //double p = Math.Atan(complex[i].Imaginary / complex[i].Real);
                //pha[i] = p * 180 / Math.PI + 90;
                coorX[i] = df * i;
            }

            return true;
        }

        private float[] Hanning(float[] input)
        {
            int N = input.Length;
            float[] hanning = new float[N];
            for (int i = 0; i < N; i++)
            {
                hanning[i] = Convert.ToSingle(input[i] * 0.5 * (1 - Math.Cos(2 * Math.PI * i / N)));
            }

            return hanning;
        }
        private void InitialPlan(float[] input)
        {
            int n = input.Length;
            // create two managed arrays, possibly misalinged
            // n*2 because we are dealing with complex numbers
            float[] fin = new float[n * 2];
            float[] fout = new float[n * 2];

            // create a managed plan as well
            mfin = new fftwf_complexarray(fin);
            mfout = new fftwf_complexarray(fout);

            // create a managed plan 
            plan = fftwf_plan.dft_1d(n, mfin, mfout, fftw_direction.Forward, fftw_flags.Estimate);
        }
    }
}
